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We consider models of identical pulse-coupled oscillators with global interactions. Previous work 
showed that under certain conditions such systems always end up in sync, but did not quantify how 
small clusters of synchronized oscillators progressively coalesce into larger ones. Using tools from 
the study of aggregation phenomena, we obtain exact results for the time-dependent distribution of 
cluster sizes as the system evolves from disorder to synchrony. 

PACS numbers: 05.45.Xt, 05.70.Ln 


In one of the first experiments on firefly synchroniza¬ 
tion, the biologists John and Elizabeth Buck captured 
hundreds of male fireflies along a tidal river near Bangkok 
and then released them at night, fifty at a time, in their 
darkened hotel room I]. As they looked on in wonder, 
they observed that “centers of synchrony began to build 
up slowly among the fireflies on the wall. In one area 
we would notice that a pair had begun to pulse in uni¬ 
son; in another part of the room a group of three would 
be flashing together, and so on.” Synchronized groups 
continued to emerge and grow, until as many as a dozen 
fireflies were blinking on and off in concert. The Bucks 
realized that the fireflies were phase shifting each other 
with their flashes, driving themselves into sync. 

Here we study stylized models of oscillators akin to 
the fireflies, in which synchrony builds up stepwise, in 
expanding clusters. By borrowing techniques used to an¬ 
alyze aggregation phenomena in polymer physics, mate¬ 
rials science, and related subjects H0, we give the first 
analytical description of how these synchronized clusters 
emerge, coalesce, and grow. We hasten to add, however, 
that the models we discuss are not even remotely real¬ 
istic descriptions of fireflies; they are merely intended as 
tractable first steps toward understanding how clusters 
evolve en route to synchrony. 

Our work is part of a broader interdisciplinary ef¬ 
fort HUH. Oscillators coupled by sudden pulses have 
been used to model sensor networks [Mmi, earthquakes 
m nn, economic booms and busts d: firing neurons 
Hi, and cardiac pacemaker cells m- Diverse forms 
of collective behavior can occur in these pulse-coupled 
systems, depending on how the oscillators are connected 
in space. Systems with local coupling often display waves 
mm or self-organized criticality Inmans. with pos¬ 
sible relevance to neural computation m and epilepsy 
m- In contrast, systems with global coupling, where ev¬ 
ery oscillator interacts equally with every other, tend to 
fall into perfect synchrony. Rigorous convergence results 
have been proven for this case [20l [22H25] . But the tech¬ 
niques used previously have not revealed much about the 


transient dynamics leading up to synchrony—the open¬ 
ing and middle game, as opposed to the end game. Ag¬ 
gregation theory offers a new set of tools to explore this 
prelude to synchrony. 

We introduce a toy model, which we call scrambler 
oscillators. It consists of N 1 identical integrate-and- 
fire oscillators coupled all to all. Each oscillator has a 
voltage-like state variable x that increases linearly ac¬ 
cording to i = 1, rising from a baseline value of 0 to 
a threshold value of 1. Whenever any oscillator reaches 
threshold, it fires and does three things, (i) It kicks ev¬ 
ery other oscillator (and every cluster of oscillators) to 
a new random voltage, independently and uniformly— 
in this sense, it scrambles the other oscillators, (ii) It 
then “absorbs” any scrambled oscillators that lie within 
a distance 1/N of threshold, by bringing them to thresh¬ 
old and thereby synchronizing with them. To avoid the 
complications that would be caused by chain reactions of 
firings, we assume that the oscillators being brought to 
threshold do not get to fire until the next time they reach 
threshold, (iii) The oscillator that fired resets to x = 0 
along with the oscillators it absorbed. 

If a cluster of j oscillators does the firing, the same 
rules apply, except that now any oscillators within a 
distance j/N of threshold get absorbed. The assumed 
proportionality to j is natural, if each member of the 
cluster contributes to the pulse strength. We study 
other plausible coupling rules in the Supplemental Mate¬ 
rial 1251.) 

The motivation for this scrambler model is that it 
leads to the simplest possible mean-field approximation. 
In the infinite -N limit, we would like clusters of every size 
to be uniformly distributed in voltage at all times. This 
convenient property would greatly ease the derivation of 
the rate equations for the cluster kinetics. As we will see 
below, the predictions that follow from this approxima¬ 
tion agree reasonably well with simulations. 

Assume the initial voltages Xi, for i = 1,..., N, are 
independent and uniformly distributed. At first, nothing 
interesting happens. The oscillators increase their volt- 
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ages without interacting. But then one oscillator reaches 
threshold and fires. The remaining oscillators get scram¬ 
bled, and perhaps some get absorbed. Then another os¬ 
cillator fires, and so on. After this process has gone on 
for a while, the system has formed synchronized clusters 
of various sizes. 

Let Nj(t) denote the number of clusters of size j at 
time t. Thus there are Ni(t) singleton oscillators, A^t) 
pairs of synchronized oscillators, Ns(t) triplets, and so 
on. The Nj are correlated random quantities. They are 
correlated because oscillators belonging to clusters of one 
size are unavailable to clusters of another size, and they 
are random because of the randomness in the initial con¬ 
ditions and the scrambling procedure. It does not seem 
feasible to understand the time-evolution of the Nj unless 
they are so large that their fluctuations from one random 
realization to another are negligible. 

So assume from now on that Nj 1 for all j and 
replace these random quantities by their ensemble av¬ 
erages. Let Cj = N~ 1 (Nj) denote the average clus¬ 
ter densities. One hopes that relative fluctuations are 
small; more precisely, N~ 1 Nj = Cj + 0(N~ 1 / 2 ). An 
even stronger assumption is that the densities of dif¬ 
ferent sub-populations are asymptotically uncorrelated: 
N~ 2 NNj = c iCj + o(N~ 1 ^ 2 ). 

These Cj allow us to define a natural disorder param¬ 
eter, given by the total density c(t) = Ej c j (*)• It mea¬ 
sures the extent of the system’s fragmentation. To see 
this, note that at t = 0 each oscillator is alone; only clus¬ 
ters of size 1 exist. Accordingly ci(0) = 1 and all other 
Cj(0) = 0 for j > 1. Hence c(0) = 1, correctly indicat¬ 
ing that the system starts out maximally fragmented. At 
the opposite extreme, as t —> oo only one giant cluster of 
synchronized oscillators exists. The system is then mini¬ 
mally fragmented: c(t) = 1/N —> 0 as N —> oo. 

To derive a rate equation for the decline of c(t ), let Ri 
be the rate at which clusters of size i fire, for i = 1,..., TV, 
and let Li be the number of clusters lost to absorption 
in each such firing. Then c = — E; RiLi. 

To find Li, recall that when a cluster of size i fires, all 
the other clusters get assigned a new voltage uniformly at 
random. Moreover, any clusters assigned to the interval 
[1 — i/N, 1) get brought to threshold and absorbed. Since 
the voltages of these other clusters are uniformly dis¬ 
tributed on [0,1], a fraction i/N of them will be absorbed. 
There arc Ej Nj clusters in total. Hence the number ab¬ 
sorbed is Li = (i/N) Ej Nj = j E c j = * c - 

The rate Ri takes more work to calculate. Since some 
clusters get absorbed, not every cluster gets the chance to 
fire. We must account for this depletion when calculating 
Ri. First consider the background rate of firing of clusters 
of size i in the absence of absorptions. In other words, 
pretend for a moment that when an j-cluster fires, it sim¬ 
ply scrambles every other cluster and restarts its own cy¬ 
cle without absorbing anyone. Call this background rate 
R Since all oscillators move with velocity u, = x t = 1, 


and since the cluster density is Cj, the corresponding 
background rate of firing is R/ = CiVi = Ci. Next, to find 
the actual Ri , we must subtract from R® the rate at which 
clusters of size i are being absorbed and hence deprived 
of their chance at firing. Call this absorption rate R 
Clusters of size i are absorbed when clusters of size j fire, 
for j = 1 ,..., N, taking a fraction j/N of the uniformly 
distributed j-clusters along with them. Since there are 
Nj clusters of size i and the j-clusters fire at rate Rj , the 
total rate at which j-clusters are being absorbed is given 

by Z>7 = Ej U/N)NjRj = Ej JCiRj = c E jRj■ 

Putting all this together gives Ri = Rj—Rj = 
Ci - Ci Ej jRj = c.i( 1 - Ej j R j)- Let fi = 1 - Ej j R j- 
Note that /3 is the same for all i, which enables it to be 
determined self-consistently, as follows. From Ri = f3ci 
we obtain (3 = 1 — Ej jRj = 1 — blow invoke 

the identity Ej j°j = j(Nj/N) = 1, which expresses con¬ 
servation of oscillators. Solving for (3 then gives (3=1/2 
and therefore Ri = Ci/2. 

Next, plug the expressions derived for Ri and Li into 
the rate equation c = — E i R iLi- The result is c = 
- Ei( c i/2)(*c) = -(c/2)Ej*c» = -c/2. Recalling that 
c(0) = 1, we conclude that 

c(t) = exp(-t/2). (1) 


Figure |T] shows this result matches simulations. 



FIG. 1: Theoretical and simulated c(t) and ci(t). Solid lines 
show theoretical curves obtained analytically (see text). Data 
points show simulation results for N = 10 4 oscillators. 

How do the individual cluster densities Ci behave? To 
derive their rate equations, note that since the voltage 
space is the interval [0,1], a segment of length TV -1 con¬ 
tains on average NcxN _1 = c clusters. In fact, the prob¬ 
ability that it contains n clusters (of any sizes) is given 
by the Poisson distribution: n„ = c n e~ c /n\. This is the 
mathematical expression of the assumption that clusters 
are distributed randomly without correlations. 

With this in mind, let us solve for Ci(t), the density 
of singletons. It is the easiest Cj(t) to analyze, since it 
can only decrease. Two mechanisms cause the loss of 
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singletons: (i) a singleton fires and absorbs a cluster, and 
(ii) a cluster fires and absorbs p > 1 singletons. 

Consider mechanism (i). From the expression for Ri 
obtained above, singletons fire at a rate R\ = c\/2. 
When they fire, they absorb any cluster lying in the volt¬ 
age segment [1 — 1/N, 1). According to the Poisson dis¬ 
tribution, the probability that this segment contains one 
or more clusters is given by 1 — e~ c . Hence singletons are 
lost by mechanism (i) at a rate (ci/ 2 ) [l — e -0 ^]. 

The loss due to mechanism (ii) is handled similarly. 
Suppose p singletons lie the interval [1 — j/N, 1) at the 
moment when a cluster of size j fires, for j = 1,..., N. 
This event happens with probability e~ JCl (jci) p /p\, and 
when it does, it consumes p singletons. The j-clusters 
fire at a rate Rj = Cj/2. Hence singletons are lost by 
mechanism (ii) at a rate 
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Summing the loss rates from (i) and (ii) gives 

— - — — (2 - 

dt ~ 2 (Z 6 ’■ 


( 3 ) 


This equation has a closed-form solution in terms of ex¬ 
ponential integrals: 

ci(i) = exp(-t + Ei(-l) - Ei(-e _t ^ 2 )), (4) 


where we have used the initial condition Ci(0) = 1. Fig¬ 
ure 1 shows good agreement between the theoretical and 
numerical ci(t). 

For i > 1, the rate equation for Ci includes gain terms 
as well as loss terms. Clusters of size i > 1 can be created 
when two or more smaller clusters coalesce, or destroyed 
when they themselves coalesce with at least one other 
cluster. The loss term is a straightforward generalization 
of that for ci, and is given by (cj/ 2 ) [2 — . 

To find the gain term, imagine that a cluster of size 
k fires. The segment [1 — k/N, 1) may contain a i clus¬ 
ters of size 1, 02 clusters of size 2, etc. This event hap¬ 
pens with probability x x ... (where we 

are using the assumption that clusters of different sizes 
are independent as well as Poisson distributed). If the 
segment contains a combination of clusters such that 
k T Oi T 2 o 2 + 3 o 3 + • • • = i, then a cluster of size i 
will form. We sum over all such combinations for a clus¬ 
ter of size k firing, and then sum over all k, to get the 
rate at which clusters of size i are created: 
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Combining the loss and gain terms, and transferring 
Cje _lc into the gain term, we finally obtain 
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We see from the sum that the equations ([ 6 ]) are recursive. 
They can be solved one by one, though not analytically, 
so we resort to numerical integration. Figure [2] shows 
that the theoretical and simulated Ci agree. 



FIG. 2: Theoretical and simulated cluster densities c 2 (t) 
through cs(t). Solid lines show theoretical predictions com¬ 
puted from numerical integration of Eq. §. Data points show 
simulation results for N = 5 x 10 4 oscillators. 

Although we cannot find all the Ci(t) explicitly, we 
can get their moments M n (t ) = V. j n Cj(t). The first 
two moments are trivial: Mq = ^ ~2jCj(t ) = c(t) and 
M\(t) = = 1) from conservation of oscillators. 

The higher moments can be obtained from a generating 
function. Let G(z,t ) = Ck{t)e kz . Then the infinite 

set of differential equations (16J) transforms into 

-—^—- + G(z, t ) = ~ C W + t)> £]■ (7) 

This equation looks neat, but it is far from trivial, as the 
right-hand side involves G in a very nonlinear manner. 
Using the identity M n (t ) = ^w| z =o, we can however 
derive the following equations for the moments: 

3 

AT 2 — -M 2 

Ms = 7 -M 3 + 3 (M 2 ) 2 ( 8 ) 

M 4 = — M 4 + 16M 2 M 3 -f- — ( 2 V/ 2 ) ^ - 

Like the Ci equations these moment equations are 
recursive and can be solved in succession. We find 

M 2 (t) = e 3t/2 

M 3 (t) = 7e n/2 - 6 e 3t (9) 

MAt) = — e 5t - 128e 9 */ 2 + ^e 15 */ 4 - 4e 27t / 8 . 

5 5 

Figure [3] plots theoretical and simulated values of the 
Mj. The agreement is good for M 2 but worse for M 3 and 
A/ 4 . This is to be expected. Each a(t) is a stochastic 
process, subject to fluctuations dominated by the chance 
formation of big clusters. Since M n (t) = j n Cj{t), the 
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FIG. 3: (Color online) Log plot of the first three nontrivial 
moments M 2 , M 3 , M 4 . Black curves, theoretical predictions 
obtained from |8|; red dots, average simulation results for 
100 realizations of N = 10 4 oscillators. 

higher moments amplify these fluctuations more and are 
therefore noisier themselves. 

Variants of the scrambler model can also be solved. 
For example, suppose that when a cluster of size j fires, 
it absorbs all oscillators within a distance kj /N of thresh¬ 
old, where k > 0 is a tunable coupling strength. Or sup¬ 
pose that the pulse strength is independent of the size 
j of the firing cluster. Both these cases are discussed in 
|26) . Nothing qualitatively new happens, but there are 
some interesting quantitative differences, e.g., in the lat¬ 
ter case the disorder parameter c(t) decays algebraically 
rather than exponentially. 

A more substantial modification would be to make the 
model deterministic. Scrambler oscillators respond ran¬ 
domly, but real oscillators respond deterministically, ac¬ 
cording to a phase-response curve m- This determinis¬ 
tic character is incorporated in Refs. [TB] [22] by assuming 
that when a cluster of size j fires, it adds a voltage pulse 
je to every other oscillator, or pulls it up to threshold, 
whichever is less. For the case where e = 1/N we show 
in [JjSJ that this deterministic system can also be well ap¬ 
proximated by the aggregation methods used above. The 
main new feature is that c(t ) and the cluster densities 
Ci(t) become piecewise linear. But their overall shapes 
still resemble those seen in the scrambler model. 

There are many avenues to explore in future work. For 
example, throughout our analysis each oscillator obeyed 
Xi = 1. Such linear sawtooth waveforms are reasonable 
for the oscillators used in sensor networks [8], but not for 
neurons or cardiac pacemaker cells. Changing the con¬ 
cavity or shape of the waveform would make the analysis 
more difficult, and might even affect the cluster kinetics, 
because clusters would no longer be uniformly distributed 
as we have assumed throughout. Could different types of 
pulses be engineered to preserve a uniform distribution? 
Or perhaps different pulses would distribute the clusters 
in a nonuniform but still tractable manner. It would 
also be interesting to study cluster kinetics in oscillator 


systems with local coupling, network structure, hetero¬ 
geneity, delays, and other realistic features. 
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